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Abstract 
A high-resolution dynamic dust source has been developed in the NASA Unified- 
Weather Research and Forecasting (NU-WRF) model to improve the existing coarse 
static dust source. In the new dust source map, topographic depression is in 1-km 
resolution and surface bareness is derived using the Normalized Difference Vegetation 
Index (NDVI) data from Moderate Resolution Imaging Spectroradiometer (MODIS). The 
new dust source better resolves the complex topographic distribution over the Western 
United States where its magnitude is higher than the existing, coarser resolution static 


source. A case study is conducted with an extreme dust storm that occurred in Phoenix, 
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Arizona in 02-03 UTC July 6, 2011. The NU-WRF model with the new high-resolution 
dynamic dust source is able to successfully capture the dust storm, which was not 
achieved with the old source identification. However the case study also reveals several 
challenges in reproducing the time evolution of the short-lived, extreme dust storm 


events. 
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1. Introduction 

Dust is one of the most abundant aerosol types in the atmosphere, playing an 
important role in the Earth’s radiation budget, cloud formation, atmospheric dynamics, 
and ocean biogeochemistry in various spatial and temporal scales (Husar et al., 2001; 
Haywood et al., 2005; Jickells et al., 2005; Forster et al., 2007; Evan et al., 2008). Mineral 
dust is also a major air pollutant that causes premature deaths by cardiopulmonary 
disease and lung cancer for the countries around the source region (Giannadaki et al., 
2014; Sprigg et al., 2014; Morman and Plumlee, 2014). The impact of dust is not limited 
to source areas but extends to larger regional or even global scales (Carlson and Prospero, 
1972; Prospero and Lamb, 2003; Kaufman et al., 2005; Chin et al., 2007; Shao et al., 
2010; Yu et al., 2012). 

The majority of global dust loading is concentrated near the sources in the 
permanent desert regions (so-called desert-belt), including North Africa, Middle East, 


and East Asia (Prospero et al., 2002; Chin et al., 2009; Huneeus et al., 2011; Ginoux et 
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al., 2012). However, dust is also emitted from semi-arid regions such as the Sahel and 
inner Mongolia, as well as from agricultural areas. Although dust aerosol generated from 
semi-arid and agricultural areas is much less than that from the major deserts, its 
importance for air quality and human health is greater at local- and regional-scales due to 
their proximity to populated areas. Correctly identifying the dust source locations and 
representing the dust storm events in numerical models are keys to estimate the impacts 
of dust on the environment and society. 

We present here the dust simulation with the NASA Unified-Weather Research 
and Forecast (NU-WRF) modeling system (Peters-Lidard et al., 2015). The objective of 
this paper is two-fold. The first goal is to describe a new, high spatial resolution (1-km) 
dynamic dust source (Sgynamic) for NU-WRF that represents an improvement of the 
existing static dust source (Static) at 0.25°x0.25° resolution (described below) currently 
available in the community WRF-Chem model. The second goal is to evaluate the NU- 
WRE model simulation of an extreme dust storm case which occurred in Phoenix, 
Arizona at 02-03 UTC July 6 (or 19-20 MST, July 5), 2011. While systematic 
observation for a severe dust storm is rare, we revisit the Phoenix dust storm which has 
been relatively better documented by observations from various platforms, including 
visual, surface radar, and surface stations (e.g., Raman et al., 2014; Vukovic et al., 2014). 
They also provide the meteorological background about the extreme dust storm. Through 
qualitative and quantitative comparisons with these direct and indirect observations, we 
discuss details about the simulated dust storm, meteorological conditions, dust source, 


and surface- and columnar intensity of the dust storm. 
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In section 2, the high-resolution dynamic dust source in the NU-WRF/GOCART 
dust emission parameterization and the model experiment setup are described. The case 
study of the Phoenix dust storm is presented in section 3. In section 4, we discuss the 


challenges in dust simulation, followed by the summary in section 5. 


2. Method 
2.1. Dust emission parameterization and source function 

The dust emission module in NU-WRF is based on the mechanisms from the 
Goddard Chemistry Aerosol Radiation and Transport (GOCART) model (Ginoux et al., 
2001). Dust emission in GOCART, assuming that the soil mobilization is proportional to 
the horizontal wind speed at near surface, is parameterized with the 10-m wind speed, the 
threshold velocity of wind erosion, and the surface condition for each dust size group 
from 0.1 to 10 um in radius (Ginoux et al., 2001, 2004; Chin et al., 2009). For each size 


group with effective radius r, dust emission flux F (ug ms’) is expressed as: 


F(r) = CS s(r) wjom(Uiom — UAT,W)); if Ujom > Uy (1) 


where C is a dimensional factor (0.4 ug s’ m~ for the current study), S is the dust source 
function or probability of dust uplifting with a value between 0 and 1, s(r) is the fraction 
of size group r within the soil, ujom is the 10m wind speed (m s'), and u; is the threshold 
velocity of wind erosion as a function of dust density, particle diameter, and surface 
wetness to account for the bonding effect between water and particles (Ginoux et al., 
2001, 2004). There are five mass size classes in the GOCART scheme with the respective 


size ranges of 0.1-1um, 1-1.8um, 1.8-3m, 3-6um, and 6-10um. The first group is clay 
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that accounts for 0.1 of the total dust mass. The balanced mass is evenly distributed to the 
remaining 4 dust size groups that are all silt. In the optical property calculations, the clay 
group is further split into four groups (0.1-0.18, 0.18-0.3, 0.3-0.6, 0.6-1) with mass 
fractions of 0.9, 8.1, 23.4, and 67.6%, respectively (Tegen and Fung, 1994). 

The topographic depression (#) and surface bareness (B) are two key parameters 
used in the GOCART scheme to calculate the dust source function ($), while other 
parameters such as soil temperature, surface wetness, and snow cover are also included in 
S calculation. The dimensionless topographic depression term H is defined as equation 2. 
H represents the probability of accumulated sediments, based on the consideration that 
dust sediments from surface erosion are accumulated in valleys and surface depressions 


(Ginoux et al., 2001; Prospero et al., 2002): 


H = (ee ‘i (2) 


Zmax—2Zmin 


where z is the altitude of a grid cell, and Za, and Zmin are the maximum and minimum 
elevations of topography in the surrounding 10°Xx10° search area. The fifth order power is 
applied to increase the topographic contrast. 

In the community WRF-Chem/GOCART, A (0.25°x0.25°) is generated with the 
GOCART scheme based on the topography and land mask from the Geophysical Fluid 
Dynamics Laboratory C360 High Resolution Atmospheric Model (~0.24° resolution) 
which are derived from the 5 min NAVY data (Ginoux et al., 2001; Putman and Lin, 
2007). The bare soil surface in the community WRF-Chem/GOCART was determined 


based on the 8 km land-cover data from the Advanced Very High Resolution Radiometer 
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(AVHRR) satellite (DeFries et al., 1998) and it does not resolve the temporal variations 
of vegetation cover. 

Recently, Kim et al. (2013) have described a method of constructing a global 
dynamic surface bareness (B) in 1°x1° spatial resolution using the 8-km spatial resolution 
AVHRR Normalized Difference Vegetation Index data (NDVI). Calculated from the 
visible (VIS) and near-infrared (NIR) radiation, NDVI reflects the state of vegetation 


over surface (Tucker, 1979): 
NDVI = (NIR-VIS)/(NIR+VIS) (3) 


MODIS NDVI has been applied for recent dust simulation studies either as a 
source masking (Vukovic et al., 2014) or as a surface vegetation fraction which is an 
input parameter for surface roughness estimation (Xi and Sokolik, 2015). In the present 
study MODIS NDVILis used to derive surface bareness following Kim et al. (2013). The 
surface is considered erodible when NDVI is below the threshold NDVI value (.e., 
NDVI). The NDVInr has been set to 0.15 taking the fact that the typical NDVI values 
are 0.05~0.10 over bare ground and the values gets larger than 0.2 during growing season 
over semi-arid region such as grass or shrub land (Huete, 1999; Zeng et al., 2000; Miller 
et al., 2006; Kim et al. 2013), such that the surface bareness B is determined as 


7 (. NDVI < NDV inp a 


0, otherwise 
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Keeping the principles of the original dynamic dust emission parameterization, 
the present study has made two major improvements. First, the degree of topographic 
depression (7) has been calculated using the U.S. Geological Survey (USGS) global 
topography map in 30 arc-second (~1-km) resolution (GTOPO30; USGS, 1996) within a 
larger search area (12°x12°). Second, the surface bareness (B) is constructed using daily 
MODIS NDVI data in 0.01° (~1-km) resolution over North America (Case et al., 2014). 
The high resolution topographic and source function better resolves the complex 
geographical variability especially over the western United States (Figure 1a and 1b). The 
MODIS NDVI for July 2011 shows a strong spatial variation ranging from 0.1 to 0.8 
(Figure lc). The erodible bare-ground (1.e., NDVI <0.15) appears over the western 


United States (Figure 1d). 


2.2. Model description and experimental setup 

The high-resolution dynamic dust source function has been implemented to the 
NU-WREF modeling system developed at NASA with collaborations with other agencies 
and institutions. NU-WREF is an observation-driven integrated modeling system that 
represents aerosol, cloud, precipitation and land processes at satellite-resolved scales at 1 
to a few km (Peters-Lidard et al., 2015). NU-WRF is built upon the Advanced Research 
WRE (ARW;; Skamarock et al. 2008) dynamical core model and the WRF-Chem (Grell et 
al., 2005), with additional NASA components that include the Goddard Land Information 
System (LIS; Kumar et al. 2006; Peters-Lidard et al. 2007), the GOCART aerosol 
modules in WRF-Chem (Tao et al., 2013), the Goddard radiation and cloud microphysics 


schemes (Tao et al. 2003; Lang et al. 2007, 2011; Shi et al. 2014), and the Goddard 
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Satellite Data Simulator Unit (G-SDSU; Matsui et al., 2013, 2014). 

For the case study, we have chosen a dust storm event that occurred in Phoenix, 
Arizona on 02-03 UTC July 6, 2011. The U.S. National Weather Service has reported 
that the Phoenix dust storm of July 5, 2011 is one of the most extreme storms in the last 
30 years (http://www. wrh.noaa.gov/psr/pns/2011/July/DustStorm.php), with an estimated 
dust front size of 160 km, depth of 18 km, and top height of 1500-1800 m. Physical 
characteristics of the dust storm and its meteorological system were analyzed by previous 
studies (Raman et al., 2014; Sprigg et al., 2014; Vukovic et al., 2014). 

In our study, a single domain encompasses 500x500 km? with the 1-km horizontal 
resolution centered at the Phoenix KIWA Weather Surveillance Radar — 1988 Doppler 
(WSR-88D) station (111.6W, 33.3N). A terrain-following, pressure-based vertical 
coordinate is used (Skamarock et al, 2008) and consists of 31 layers with a model top 
pressure of 50 hPa. The model integration time step is set to 3 seconds. The model was 
simulated for 48 hours, from 00 UTC 5 July, 2011 to 00 UTC 7 July, 2011. Table 1 
summarizes the NU-WREF configuration options selected for various atmospheric 
processes. Initial and lateral boundary conditions for meteorological variables were 
obtained from the NCEP Global Forecasting System (GFS). The NU-WRF/GOCART 
simulations were also conducted for anthropogenic aerosols using GOCART aerosol 
model and its contribution to PM10 is less than a few wg m° during the storm. For the 
case study, the empirical dimensional factor (C) is set to 0.4, emitted dust is equally 
distributed to the lowest 5 model layers which extend to about 500 m above ground level 


(AGL), and the cutoff soil moisture factor (gyet) 1s set to 0.35. In the present case study, 
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we only consider dust aerosol because of the negligible level of other aerosols during the 


dust storm. 


3. Results of the case study 
3.1. Phoenix dust storm on July 5, 2011 

The U.S. National Weather Service has reported that strong thunderstorms 
developed east of Tucson, AZ during the afternoon local hours of July 5, 2011 
(http://www.wrh.noaa.gov/psr/pns/2011/July/DustStorm.php). The storms intensified as 
they progressed west into the Tucson Metropolitan Area, producing downburst winds in 
excess of 30 ms”. The leading edge of these strong outflow winds moved to the 
northwest at 45 to 65 km hr”. The first dust wall in Casa Grande (which is located at 
about 75 km southwest of Phoenix) was reported to NWS Phoenix by 06:30 PM MST 
July 5 (01:30 UTC July 6). The unstable atmosphere with a convective available potential 
energy of 786 J kg! marked at the Tucson weather station at 12 UTC July 5, 2011 
suggests conditions are favorable for a strong thunderstorm. 

The thunderstorm near Tucson that initiated the Phoenix dust storm is associated 
with the North American Monsoon (NAM), when the synoptic scale wind and rainfall 
shifts in the summer over Mexico and the southwestern U.S. (Douglas et al., 1993; 
Adams and Comrie, 1997). Briefly speaking, the NAM circulation pattern typically 
develops in late May or early June over southwest Mexico, moving to northwest Mexico 
in mid to late June, and to the southern U.S. in early July 
(http://www.wrh.noaa.gov/twc/monsoon/monsoon_NA.php). In the NAM, the low level 


moisture is transported primarily from the Gulf of California and eastern Pacific. The 
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upper level moisture is transported mainly from the Gulf of Mexico by easterly winds. 
Since the lower level moisture flow is not as persistent, the state of the upper level 
atmosphere is important for thunderstorm development on a given day. In addition, the 
southwestern U.S. was experiencing a moderate to extreme drought in 2011 according to 
the U.S. National Drought Monitor (http://www.drought.gov/drought/). The barren 
surface resulting from the drought provides more favorable conditions for dust emission. 

The Phoenix KIWA radar (111.6°W, 33.3°N, 412 m Mean Sea Level, MSL) was 
operating during the Phoenix dust storm event. Although the main application of Next 
Generation Weather Radar Level-III (NEXRAD Level-II]) is for weather analysis and 
forecasts, previous studies (e.g., Raman et al., 2014; Vukovic et al., 2014) have used it to 
probe the location, area, and motion of dust storm since more reliable and continuous 
dust storm observations are absent. The radar data is available with 5-minute frequency at 
the National Centers for Environmental Information NEXRAD online data inventory 
(http://www.ncdc.noaa.gov/nexradinv/index.jsp). In Figure 2, we show hourly radar 
reflectivity, co-polar correlation coefficient, and base velocity from the NEXRAD Level- 
II at 0.5° elevation angle between 01:30 to 03:32 UTC. The radar reflectivity is the 
returned signal within a sampling volume to the radar station. The co-polar correlation 
coefficient is the correlation between the backscattered horizontal and vertical polarized 
signals ranging zero (i.e., non-spherical shape) and one (i.e., spherical shape). The base 
velocity is the measure of the radial component of the wind either negative (1.e., toward) 
or positive (i.e., away) values from the radar. Although limited, the strong dust storm can 
be identified with combined analysis of the radar and surface observations. 


At 01:30 UTC July 6, the radar showed well-defined bow shape in the 


10 


229 


230 


Zol 


Zoe 


230 


234 


235 


236 


257 


238 


Zo9 


240 


241 


242 


243 


244 


245 


246 


247 


248 


249 


250 


251 


southeastern area of the KIWA as marked with an arrow in Figure 2a. The bow shape is 
characterized with weak reflectivity (<20 dBZ), low co-polar correlation coefficient 
(<0.8), and a stronger wind (20 knots or 10 m s) toward KIWA. Combined with the 
visual observations by NWS reports, previous studies (e.g., Raman et al., 2014; Vukovic 
et al., 2014) have suggested that the bow shape is the dust storm that hit the Phoenix 
Metropolitan area. At 02:31 UTC, the dust front moves toward the KIWA radar station 
with the extended front size. The wind front passed the radar station and the sign of the 
base radial velocity is now positive. At 03:32 UTC, the dust storm continues moving 
toward the northwest and the front is located inside of the Phoenix Metropolitan area 
(i.e., the area of the white rectangle). Key information from the radar analysis includes: 
(1) the maximum of the base velocity is larger than 23 ms” (or 45 knots) and its origin is 
located at the southeast of Phoenix near Tucson; (2) the front moves toward northwest 
with a speed in 50-60 km hour’; and (3) the major dust storm area covers Phoenix at 02- 


03 UTC and its front expands more than 150 km. 


3.2. NU-WRE simulation of the Phoenix dust storm 

From this section, our case study domain covers 500 km x 500 km centered on the 
KIWA radar station. The 1-km topography map shows that the northeastern region of the 
domain mainly consists of mountains higher than 1500 m, while the southwestern region 
consists of lower terrain with heights below 500 m (Figure 3a). H is low over the 
mountain regions (<0.1) but it is higher over the coastal regions and basins (>0.2) due to 
the inverse relationship with topography (Figure 3b). The surface conditions are 


characterized by the arid southwestern region with the NDVI below 0.15, and the 
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vegetated northeastern mountainous region with the NDVI larger than 0.2 (Figure 3c and 
3d). The resulting high-resolution dynamic dust source function (Sgynamic) Covers most of 
the southwestern basin region ranging from 0.05 to 0.3, but it excludes most of the 
northeastern mountains (Figure 3e). In contrast, the static dust source (Ssiaric) does not 
show the detailed geographical structure and its values (mostly below 0.05) are much 
smaller than those in Sgynamic for the same geographic areas (Figure 3f). 

The horizontal 10-m wind field (W10m) from NU-WRE is plotted every 2 hours 
during the dust storm (Figure 4). A strong wind area begins to form at 21 UTC July 5, 
2011 and it develops a clear wind gust front two hours later (23 UTC) over the southern 
region of the domain. The magnitude of the wind gust continues to intensify and the 
maximum wind speed rapidly moves toward northwest (01 UTC July 6). The area with 
strong wind passes through the Phoenix Metropolitan area at 03 UTC. Then it continues 
to move towards the northwest until it weakens in the next four hours. The maximum 
simulated wind speed in the study domain during the event is larger than 20 m s'. The 
model captures the initial location and fast motion of the storm observed by the radar 
shown in the previous section. On the other hand, the simulation also shows that the 
maximum strength of the storm is located further west than the radar observation (i.e., 
Figure 2). 

The temperature contours and y-z component of the cross section wind vector that 
passes through Phoenix (i.e., shown in Figure 3a) are plotted during the dust storm period 
from 21 UTC July 5 with a 2-hour interval (Figure 5). At 21 UTC July 5, the vertical 
atmospheric structure is characterized with the gradual decreasing change of temperature 


from near surface (> 30°C) to upper air (15°C at 3 km MSL). The north-south wind 
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vector component is not strong at the time, but the presence of updraft near 31°N 
indicates that the thunderstorm is beginning to develop. At 23 UTC, the temperature of 
the lower atmosphere in the south of 32°N significantly drops from 33°C to 24°C as a 
result of the rain-cooled downburst from the upper atmosphere. At 01 UTC July 6, the 
strong horizontal blowouts of cold air continue to progress to north. At 03 UTC the 
strong wind extends to 34°N with the curl shaped wind pattern in the lower atmosphere 
of 1-3 km MSL. The intensive downburst is dissipated at 05 UTC (Figure 5). Our result is 
consistent with the NWS report that the explosive horizontal outflow during the Phoenix 
dust storm was initiated by downburst generated by the thunderstorm, which occurred 
near the Tucson Metropolitan Area. 

Hourly surface meteorology data over the study domain is available from 
NOAA's National Centers for Environmental Information (NCEI) 
(https://www7.ncdc.noaa.gov/CDO/cdo). We analyze the time series of meteorological 
fields of wind speed, temperature, relative humidity, and surface pressure at three airport 
stations of Tucson (110.96W, 32.13N, 777m), Casa Grande (111.77W, 32.95N, 446m), 
and Phoenix (112.00W, 33.43N, 337m) and compare them with the NU-WRF model 
(Figure 6). Located from South to North along the storm track, the station data provides 
useful insight of the surface meteorological conditions during the storm passage. The 
most noticeable result form the observation is the rapid change of meteorological fields 
with the arrival of the storm. At Tucson, for example, the observation shows that the 
rapid change occurs between 23 UTC July 5 and 01 UTC July 6. Temperature decreases 
from 36.7 °C to 21.7 °C and wind speed increases from 5.8 ms” to 8.9 ms"! with the 


arrival of the storm. The increase of relative humidity from 21% to 82% is explained with 
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the large decrease of temperature. Surface pressure does not vary much ranging from 920 
to 930 mb. Similar sudden changes appear in the Northern stations of Casa Grande (01- 
03 UTC July 6) and Phoenix (03-05 UTC July 6), but 2 to 4 hour later than Tucson. At 
Casa Grande, the temperature and wind speed are changed from 38.9 °C to 23.9 °C and 
from 1.3 ms"! to 15.2 ms”, respectively. At Phoenix, the temperature and wind speed are 
changed from 37.8 °C to 27.2 °C and from 3.6 ms’! to 8.9 ms", respectively. NU-WRF 
model captures the magnitude and pattern of the observation, showing that it can 
reproduce the storm and its evolution. However the comparison also shows that the 
simulated storm is moving faster than observation resulting 1 or 2 hour earlier storm 
arrival at Casa Grande and Phoenix. Daily accumulated precipitation was 25.6 mm at 
Tucson station, but no precipitation is reported at Phoenix station or negligible at Casa 
Grande station (<1 mm). 

The dust emission is plotted in Figure 7, and is mainly controlled by W10m since 
the dust emission in the NU-WRF/GOCART is proportional to the 3™ order of W10m 
(Eq. 1). The amount of dust emission exceeds 100 pg ie during peak dust storm hours 
of 01-03 UTC. In contrast to the original soil moisture threshold value of gwet < 0.5, a 
reduced threshold values (1.e., wer < 0.35) was used in the current simulation to achieve a 
better agreement with the radar observation. As a result dust emission over the 
southwestern region of the domain (i.e., southwest of 113°W, 33°N) is substantially 
suppressed during the dust storm period. The time-evolution of the surface dust PM10 
(dust size is less than 10 um) concentrations is quite similar to that of dust emission 
(Figure 8). The model simulated PM10 over the Phoenix area (i.e., inside the black 


square) is less than 100 pg m® most of time at 01 UTC July 6, 2011. When dust storm 
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reaches the Phoenix at 03 UTC, the PM10 drastically increases to reach 4000 ug m°, 
After the dust storm at 05-07 UTC, the PM10 concentration gradually reduces but still 
remains much higher than that before the dust storm. In section 4, we will discuss the 
apparent tardy decay of dust concentrations simulated by the model. 

The cross sections of dust concentration and the y-z component of wind vector at 
112°W during the dust storm are plotted in Figure 9. At 21-23 UTC July 5, model- 
simulated dust concentrations are relatively low because of the low emission rates (Figure 
7). At 01 UTC July 6, high dust concentration (700-1000 pg m”) first appears in the 
latitudinal zone between 32.3°N and 32.8°N that agrees with the horizontal dust surface 
concentration field in Figure 8. The model shows that at 2 km above the ground the dust 
concentration is about 100 yg m™. The high dust concentration area rapidly moves north 
following the strong horizontal wind gust. At 03 UTC, the front of dust moves about 60 
km hour? reaching to 33.9°N. At the same time, the lofted dust layer (>1 km above 
surface level) is found behind the dust storm between 32.6°N and 33.9°N. At 05 UTC, 
dust storm front continues moving north but its moving speed has reduced by about half. 
However, the surface dust concentration remains high (> 1500 pg m”) at 33-34.5°N after 


the storm at 05-07 UTC, as shown in Figure 8. 


3.3. Comparisons of model simulation with observations 

The simulated dust at surface level is compared with the air quality data from the 
U.S. Environmental Protection Agency’s Air Quality System obtained from the EPA’s 
AirData website (http://www.epa.gov/airdata/) (Figure 10). Across the state of Arizona, 


all of the 13 sites with clear dust storm signals are in and around the Phoenix 
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Metropolitan area in the range between 113.3°W, 32.0°N and 110.4°W, 34.3°N (Le., the 
black box in Figure 3a). All EPA PM10 observations indicate a sharp dust peak that 
occurred in the two-hour window between 02 to 04 UTC on July 6 with a maximum at 03 
UTC, although the PM10 magnitudes vary by location from 1946 to 6348 pg m”. The 
model captures the observed peak events at most sites with an averaged value over 13 
sites (2968 ug m”) similar to the observations (2505 ug m”) and the average correlation 
coefficient of 0.63. In contrast, a run with the static dust source has simulated only a 
negligible amount of surface PM10 concentration (81~258 j1g m™) in those EPA stations 
during the same period (i.e., Figure 10). After the dust storm, the observations show a 
rapid decrease of dust concentrations at all sites to the pre-storm levels, but the dust 
concentration in the NU-WRF model remains elevated. This after-peak high bias in the 
model could be caused by various physical reasons such as the location and progress of 
dust storm and uncertainties in dry deposition or emission processes. Further 
investigation with various model runs has found that dust emission from south of the 
Phoenix area for 2 hours from 0300 UTC July 6 is most responsible to the high bias in the 
NU-WRE model. A sensitivity simulation that turned off dust emissions after 0300 UTC 
indeed removes the high dust residual after the storm and improves the correlation 


coefficient (r=0.89) and other statistics (Figure 11). 


4, Discussion 
Although the case study shows the high-resolution dynamic dust source 
considerably improves dust modeling, it also illustrates several outstanding challenges in 


dust emission processes in the NU-WRF/GOCART model: 
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(iii) 


The curly motion of the outflowing dust front in the downburst produces the 
wall of dust, which reached higher than 1.5 km in the 2-3 hour time span. 
However, the advection/convection scheme in NU-WREF could not resolve the 
rapid vertical transport of the high dust wall, leaving most of the dust in the 
lowest levels. In the present case study, we equally distributed the emitted 
dust to the lowest 5 model layers (which are about 0-500 m above the ground) 
to better resolve the vertical distribution of the simulated dust. It is necessary 
to consider a better mechanism to represent the vertical distribution of emitted 
dust in the “haboob” events. 

Dust emission is inhibited when gwet is larger than the threshold. While gwet 
values are regionally dependent from 0.35 to 0.5 in global GOCART 
modeling studies (e.g., Kim et al., 2013; Chin et al., 2014), the threshold gwet 
of 0.35 is better compared to the observations in the present case study. More 
robust constrain on the threshold gwet is necessary in future studies. 

The dimensional factor (C) in the global-scale GOCART is set to 1 ug sm 
(Ginoux et al., 2001; Kim et al., 2013). But C values are highly case- 
dependent, and are in the range from 0.65 to 22 zg s’ m™ in previous WRF- 
Chem studies with the static dust source (Zhao et al., 2010; Bian et al., 2011; 
Alizadeh Choobari et al., 2013; Kumar et al., 2014). In the present extreme 
dust storm case study, the C value was set to 0.4 wg s* m> to achieve better 
agreement with the observations. A more generalized method of setting C 


values is required in future studies. 
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The case study also showed that simulating the correct wind field for calculating 
dust emission is very important but challenging. Although the NU-WRE used a realistic 
meteorology (i.e., meteorological fields from reanalysis or model) to initialize simulation 
and force the lateral boundaries, the location and time evolution of the wind storm within 
the regional domain are still problematic. For example, the center of the outflowing wind 
storm in our simulation is positioned too far west compared to the radar. We conducted 
20 sensitivity runs with different modeling setup and configuration options by varying 
domain related configurations (initial time, domain nesting, horizontal- and vertical- 
resolutions) and physics related schemes (planetary boundary layer schemes, longwave- 
and shortwave-radiation schemes, land-cover, land surface model, initial meteorology 
input, aerosol radiative feedback, and data assimilation). While most runs successfully 
captured general characteristics of the Arizona dust storm, no simulation successfully 
captured the exact timing, location, and evolution of the storm. 

Similarly to our study, a previous study simulated the Arizona dust storm using 
the Nonhydrostatic Mesoscale Model on E grid-Dust REgional Atmospheric Model 
(NMME-DREAM) with 4 km horizontal resolution (Vukovic et al., 2014). The model 
successfully simulated the position of the front in space and time and horizontal and 
vertical distribution of dust. Using MODIS NDVI, they have highlighted the importance 
of vegetation masking to improve dust simulations. Similarly to our modeling study, they 
also showed some challenging issues in the model results, such as the 1 hour late storm 
arrival time in Phoenix, underestimations of surface PM10 (<2500 pg m>), and a strong 


residual dust 4 hours after the storm (in their Figure 9). 
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5. Summary 

In the present study, we have developed a high-resolution dynamic dust source in 
the NU-WRF model. The source function is calculated from the 1-km topography map 
and from the surface bareness derived from the dynamic surface vegetation information 
from MODIS at 1-km resolution. A case study simulating an extreme dust storm occurred 
in Phoenix, Arizona in 02-03 UTC July 6, 2011 has demonstrated that the new high- 
resolution dynamic dust source better captures the complex topographic distribution 
pattern and it simulates the dust storm better than the previously used lower resolution, 
static dust source in this case. Although there is some discrepancy, the model captures the 
initial location and fast motion of the storm observed by the radar. 

NU-WRE surface dust PM10 is compared with the 13 station data in the EPA’s 
Air Quality System network. The time series analysis shows that NU-WRF can 
successfully simulate the progress of the Phoenix dust storm (R=0.63) and its magnitude. 
At the peak hour at Phoenix (03 UTC July 6), the PM10 drastically increased with the 
observed and simulated station means of 2968 and 2505 pg m, respectively. 
Significantly elevated dust PM10 values after the dust storm (e.g., at 07 UTC) simulated 
by NU-WREF were found to be due to excess dust emission near the Phoenix region 
between 03-04 UTC, when the actual dust storm had already passed the city. 

The NU-WREF model with the new high-resolution dynamic dust source is able to 
capture the Phoenix dust storm, which was not possible using the old static sources. 
However, the case study also has revealed several issues in the NU-WRF/GOCART 
model to reproduce the rapid change of surface concentrations during the event. The NU- 


WRE model could not exactly place the location of outflowing winds against radar, even 
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after 20 additional sensitivity runs with different configurations and physics options. This 
highlights that simulating accurate meteorology and wind fields is highly important for 


dust storm prediction but it is also a challenging task. 
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Table 1. NU-WRF physics setup 


Processes Name Setup 
Longwave radiation Goddard 
Shortwave radiation Goddard 
Surface layer Monin-Obukhov 
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Boundary layer YSU 
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Figure 1. (a) Surface elevation (m), (b) degree of topographic depression (H), (c) MODIS 
NDVI for July 2011, and (d) location of surface bareness (B) over North America. Black 
square in (a) is the model domain for the case study. 
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Figure 2. (a) Radar reflectivity in dBZ, (b) co-polar correlation coefficient, and (c) base 
velocity in knots (1 knot=0.51ms’) at 0.5° elevation angle from the NEXRAD Level-III 
during the dust storm event between 0130-0332 UTC July 6, 2011. The location of the 
KIWA radar station (111.6°W, 33.3°N, 412m) is marked as star symbol. The wind 
direction is toward the radar when the base velocity is negative or vice versa. Thick 
dashed-lines indicate the location of dust storm front. Domain covers 
[113.2°W~110.61°W; 32.0°N~34.1°N] and the white rectangle is the Phoenix 
metropolitan area which is mixed with some shrublands and grasslands. 
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Figure 3. Surface elevation (m), topographic depression (H), NDVI, bareness (B), 
dynamic source function (Spynamic), and static source function (Sstatic) over the dust storm 
case study domain on July 5, 2011. Phoenix Metropolitan is located within the black box. 
KIWA and KEM<X are radar stations at Phoenix and Tucson. 
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Figure 4. NU-WRF 10-m wind speed (ms”) and vector over the study domain from 21 
UTC July 5 to 07 UTC July 6, 2011. Black square indicates the NEXRAD radar domain 
shown in Figure 2 and white rectangle is the Phoenix Metropolitan area. Thick dashed- 
lines indicate the location of dust storm front from radar observation. 
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Figure 5. NU-WRF vertical cross section of temperature (°C) and the v-w component of 
wind vector (ms"') from 21 UTC July 5 to 07 UTC July 6, 2011. The cross section is 
along 112°W as shown in Figure 3a. The white area is the topography height. 
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(b) Casa Grande (111.77W,32.95N) 
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Figure 6. Meteorological variables from weather station measurements (dotted 
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NU-WREF model (solid-line) from 20 UTC July 5 to 09 UTC July 
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Figure 7. Same as Figure 4 except for dust emission (ug ms"). 
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Figure 8. Same as Figure 4 except for surface dust PM10 concentration (ug m°). 


w 
b 
CMI ms TD 


20 


Height (km) 


w 
Nilateteiat | 11 5 | 


Height (km) 


Ee 
i) 


5 50250 1000 2000 3000 5000 


Figure 9. Same as Figure 5 except for surface dust PM10 concentration (ug m™). 
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Figure 10. Time series of hourly mean observed PM10 and modeled dust PM10 at 13 
EPA stations from 20 UTC July 5 to 09 UTC July 6, 2011. Static source results are 
shown in blue dashed lines. The panel on the bottom-right is the mean of 13 station 


values. 


$1[112.14W; 33.48N] 


05:20 05:23 06:02 06:05 06:08 
Time (dd:hh) 


S6[111.88W; 33.30N] 


2 


$2[112.19W; 33.57N] 


5:20 05:23 06:02 06:05 06:08 
Time (dd:hh) 


$7[111.72W; 33.31N] 


05:20 05:23 06:02 06:05 06:08 


$3[112.05W; 33.46N] 


O(ez)=8078 B= 177 
M(3z)=2555 E= 744 
R= 0.68 


Time (dd:hh) 


$8[112.14W; 33.41N] 


S4[112.12W; 33.46N] 


B= 288 


Mlae)=2784 E= 464 


| R= 0.89 


05:20 05:23 06:02 06:05 06:08 
Time (dd:hh) 


$9[112.34W; 33.64N] 


S5[112.08W; 33.40N] 


20 05:23 06:02 06:05 06:08 
Time (dd:hh) 


$10[112.62W; 33.37N] 


3000 7 o(az)=2067 B= 10 6000 +o(gz)=5789 B= 180 3000 7OGz)=1046 Ba 254 3000 7-O(6z)=2757 B= 177 3000 7 O7)=2186 B 
2500 4 M(3z)=1514 5000 4 M(3z)=185 L 2500 2500 4 M(3z)=1664 
a es oc £ R= 0.83 
£ 2000 £ 4000 4 F € 2000 —£ 2000 4 
= 1500 = 3000 = 1500 = 1500 4 
oOo Oo Oo Oo 
S 1000 S 2000 4 - 5S 1000 S 1000 4 
a a a oa 
500 4000 4 - 500 500 4 
0 0 0 0 
05:20 05:23 06:02 06:05 06:08 05:20 05:23 06:02 06:05 06:08 05:20 05:23 06:02 06:05 06:08 05:20 05:23 06:02 06:05 06:08 05:20 05:23 06:02 06:05 06:08 
Time (dd:hh) Time (dd:hh) Time (dd:hh) Time (dd:hh) Time (dd:hh) 
$11[112.29W; 33.69N] $12[112.12W; 33.43N] $13[112.10W; 33.50N] Station Mean 
3000 7O(3z)=379 B= 369 3000 7 O(a7)=1974 B= 293 O(Gz)=2568 B= 206 
2500 M(3z)=2712 5000 M(3z)=2508 E= 396 
"E 2000 E 
2 4500 2 —OBS 
Oo Oo 
S 1000 S 
a a SS = 
a MOD 
0 


05:20 05:23 06:02 06:05 06:08 
Time (dd:hh) 


05:20 05:23 06:02 06:05 06:08 
Time (dd:hh) 


05:20 05:23 06:02 06:05 06:08 
Time (dd:hh) 


05:20 05:23 06:02 06:05 06:08 
Time (dd:hh) 


Figure 11. Time series of hourly mean observed PM10 and modeled dust PM10 at 13 


EPA stations from 20 UTC July 5 to 09 UTC July 6, 2011. The panel on the bottom-right 
is the mean of 13 station values. In the model, emission is prohibited for 0300-0500 UTC 
July 6, 2011. 


Highlights: 
e A high-resolution dynamic dust source has been developed. 
e New dust source better resolves the complex topographic distribution. 


e Acase study is successfully conducted with a strong dust storm in NU-WRF. 


